#ifndef LAMINAR_FLOW_H
#define LAMINAR_FLOW_H

#include <vtkDataSetReader.h>
#include <iostream>
#include <vector>

namespace cfd {
using namespace std;

//typedef vector<double> Vec;
//typedef vector<Vec> Mat;

template<typename T>
class LaminarFlow {
  /**
   * _it         iteration
   * _n          number of annuli
   * _mu_const   dynamic viscosity
   * _L          length of tube
   * _p_0        pressure at z = 0 (input)
   * _p_L        pressure at z = L (output)
   * _r          vector of radii (from the center to the outside)
   * _grad_p     gradient of pressure <math>-\frac{p_l-p_0}{L}=\frac{dp}{dz}</math>
   * 
   **/

  /**
   * Matrix A is tridiagonal and saved in the form of a vector of a vector with 3 components. 
   *   ( x x x              )
   *   (   x x x            )
   *   (     x x x          )
   *   (       ...          )
   *   (            x x x   )
   *   (              x x x )
   * the first and last column are "ghost rows" without any meaning.
   *
   **/

  private:
    const int _n;
    int _it;
    const T _mu_const, _L, _p_0, _p_L, _c;
    T _grad_p;

    vector< vector<T> > _A;
    vector<T> _w, _b;
    
    vector<T> _mu, _h, _r;   

    void resolveLES();
    void printA();
    void showVector(vector<T> *vec);

  public:
    // 1 Initialization
    LaminarFlow(int n=100, T mu_const=1.0, T L=1.0, T p_0=1.0, T p_L=0.0, T c=10.0);

    ~LaminarFlow();
    
    // 2 One step
    void oneStep();
    void plot();
};

} // namespace cfd

#endif // LAMINAR_FLOW_H

